require(dhlvm)

issue.labels <- c("消費税率10%超", "年金医療費抑制", "産業保護抑制", 
                  "PB均衡先送り", "量的金融緩和続行", 
                  "防衛力強化", "敵基地攻撃", "対北朝鮮圧力", 
                  "非核三原則", "辺野古移設", "小さな政府", 
                  "公共事業", "財政出動", "消費税率下げ", 
                  "累進課税強化", "法人税率上げ", "外国人受け入れ", 
                  "治安優先", "処理水放出", "夫婦別姓", 
                  "同性婚", "理解増進法", "日米同盟強化", 
                  "対中認識", "格差より経済", "国内産業保護", 
                  "原発廃止", "財政赤字", "合区解消", "憲法改正", 
                  "自衛隊", "集団的自衛権", "環境権", 
                  "プライバシー権", "知る権利", "家族条項", 
                  "教育の充実", "都道府県代表", "首相公選制", 
                  "財政の健全性", "地方分権", "憲法裁判所", 
                  "改正発議要件緩和", "緊急事態条項", "解散権制約")

#### データの準備 ####
## 東京大学谷口研究室・朝日新聞2021年衆院選政治家調査データ
elite.data <- read.csv("2021UTASP20211126_excludedQ11.csv")
# 全項目無回答者を除外
elite.data <- subset(elite.data, RESPONSE == 1)
# # 争点態度のデータを抽出
# issue.data <- elite.data[, 51:80]
# # 欠損値を6に置換
# issue.data[elite.issue.data == 99] <- 6

# 争点態度のデータを抽出
issue.data <- elite.data[, 51:95]
# 欠損値を6に置換
issue.data[issue.data == 99] <- 6
# 改憲項目に関しては欠損値を3に置換
for (i in 1:15) {
  issue.data[, 30 + i] <- ifelse(issue.data[, 30 + i] == 66,
                                 ifelse(issue.data[, 30] == 6, 3, 1),
                                 ifelse(issue.data[, 30 + i] == 6,
                                        3, issue.data[, 30 + i] + 1))
}

#### 予備的分析のためのLDA-Sのパラメータの推定 ####
y <- as.matrix(issue.data)
G <- max(elite.data$PARTY)

## K = 2, ..., 10としてパラメータを推定
preliminary.result <- list()
for (i in 2:10) {
  set.seed(i)
  K <- i
  eta <- list()
  for (j in 1:30) {
    eta[[j]] <- matrix(1, K, 6)
  }
  for (j in 31:45) {
    eta[[j]] <- matrix(1, K, 3)
  }
  alpha <- matrix(1, G, K)
  preliminary.result <- hlcModel(y, elite.data$PARTY, eta, alpha, 10000, 2000, 1)
  save(preliminary.result, file = paste0("preliminary_result_", i, ".Rdata"))
}

# BICの比較
BIC <- rep(NA, 9)
for (i in 2:10) {
  load(paste0("preliminary_result_", i, ".Rdata"))
  posterior <- posteriorMeans(preliminary.result)
  BIC[i - 1] <- bic(as.matrix(issue.data), elite.data$PARTY,
                    posterior$pi, posterior$beta)
}
which.min(BIC) + 1
# 7グループモデルを採用

#### 政策クラスターを分類する上での争点の重要度の評価 ####
load("preliminary_result_7.Rdata")
# パラメータの事後平均
preliminary.posterior <- posteriorMeans(preliminary.result)
names(preliminary.posterior$beta) <- issue.labels

## 付録C.2のD_k（式(5)）の計算
# D_kを計算する関数
dist.fun <- function(x, member) {
  n <- nrow(x)
  dist <- matrix(NA, n, n)
  weight <- matrix(NA, n, n)
  for (i in 1:n) {
    for (j in 1:(n - 1)) {
      if (i == j) {
        break
      }
      dist[i, j] <- norm(matrix(x[i, ] - x[j, ]))
      weight[i, j] <- sum(member[c(i, j)]) / sum(member)
    }
  }
  weighted.mean(dist, weight, na.rm = TRUE)
}

# D_kの計算結果（表A2）
item.distance <- sapply(preliminary.posterior$beta[1:30], dist.fun, 
                        member = table(preliminary.posterior$z_assign))
round(rev(sort(item.distance)), 3)